Constraint Damping in First-Order Evolution Systems for Numerical Relativity 
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A new constraint suppressing formulation of the Einstein evolution equations is presented, gen- 
eralizing the five-parameter first-order system due to Kidder, Scheel and Teukolsky (KST). The 
auxiliary fields, introduced to make the KST system first-order, are given modified evolution equa- 
tions designed to drive constraint violations toward zero. The algebraic structure of the new system 
is investigated, showing that the modifications preserve the hyperbolicity of the fundamental and 
constraint evolution equations. The evolution of the constraints for pertubations of flat spacetime 
is completely analyzed, and all finite-wavelength constraint modes are shown to decay exponentially 
when certain adjustable parameters satisfy appropriate inequalities. Numerical simulations of a 
single Schwarzschild black hole are presented, demonstrating the effectiveness of the new constraint- 
damping modifications. 
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I. INTRODUCTION 

Numerical relativity has recently undergone a revolu- 
tion. Multiple research groups, using a variety of math- 
ematical and computational formalisms, have produced 
consistent pictures of the late mspiral and coalescence 
of binary black hole systems 0, 0, i, i, i, 1, 0, H, a 
goal that until recently seemed remote. The community 
is opening a theoretical window on issues fundamental 
to gravitational wave astrophysics, but much work still 
remains to be done. 

State of the art simulations have provided gravitational 
waveforms due to the last several orbits of a binary black 
hole system, and its eventual coalescence and ringdown. 
This is an extraordinary achievement, but numerical rel- 
ativity may need to handle well over a dozen orbits ac- 
curately before a seamless transition can be made from 
post-Newtonian analysis. This requires an exceptionally 
stable evolution scheme, and also an extremely efficient 
one, if the simulations are to be done in a timely manner. 

Numerical simulations can become unstable for a va- 
riety of reasons, some purely numerical (such as a poor 
choice of algorithm) , some purely mathematical (such as 
ill-posedness of the continuum mathematical problem), 
and some a combination of the two. The subject of this 
paper is an instability of this last type: the exponential 
growth of the constraint fields under free evolution - an 
instability of the continuum evolution equations, seeded 
by numerical errors. 

A variety of methods exist to deal with such insta- 
bilities. One very well established method is known as 
constrained evolution @, 01, HO, 03, EE El EE EH, in 
which some subset of the dynamical fields are integrated 
in time using the evolution equations, and others are ob- 
tained by solving the constraints after each time step. 
This separation of the fields into an "evolved" family and 
a "constrained" family is normally guided by some sort 
of symmetry. A related method, known as constraint 
projection [IE EH , places all fields on an equal footing, 
freely integrating everything using the evolution equa- 



tions, then periodically "projecting" the fields down to 
the constraint-satisfying subset of the solution space. It 
appears that this method can be quite robust in practice, 
but it can also be technically demanding, requiring the 
repeated solution of nonlinear elliptic equations. 

A preferable approach for dealing with these instabili- 
ties, whenever possible, is to remove them from the evo- 
lution system at the continuum level, before they reach 
the numerical code. This type of effort is often referred 
to as "constraint damping." It is possible to change the 
evolution equations without changing the physics they 
represent. For instance, coordinates can be chosen freely 
on the simulated spacetime, and indeed, careful choices 
of gauge have been shown to have a strong effect on the 
stability of simulations, particularly in the BSSN sys- 
tem [HI EE El E3- Another, perhaps more drastic 
method to stabilize constraints involves extending the 
family of evolved fields. It is possible, with some care, 
to introduce fields whose evolution will naturally lead to 
the presence of friction terms in the implied constraint 
evolution system. Systems of this form have come to 
be referred to as "A systems" in the relativity literature, 
after the pioneering work of Brodbeck et al. [IE E3| • 

The constraint damping of Ref. [23| utilizes the free- 
dom to substitute constraint equations into the evolution 
equations. In the physical situation where the constraints 
are all satisfied, the equations are unchanged. However 
in the numerical situation where the constraint satisfac- 
tion is at best approximate, these substitutions can have 
a profound effect on the structure of the evolution sys - 
tem and the stability of its constraints. In Ref. [25j |. 
Kidder, Scheel and Teukolsky added terms, proportional 
to the constraints, to a first-order representation of the 
Arnowitt-Deser-Misner (ADM) evolution equations [26j 
(starting from the form advocated by York [27j )■ With 
these modifications, they were able to change the princi- 
pal part of the evolution system. When the free parame- 
ters satisfy certain inequalities, the evolution system be- 
comes strongly, or symmetric, hyperbolic [25l. l28l. l2fj|. [3f| . 
The purpose of the present paper is to generalize the KST 
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systems even further, to add more terms proportional to 
constraints into the evolution equations, but now with 
the goal of damping the constraints while preserving hy- 
perbolicity. It will be shown that this goal can largely be 
achieved, stabilizing all the constraints of the KST sys- 
tems, without the need to introduce extra fields as in the 
"A system" approach. 

The possibility of constraint damping along these lines 
is now widely seen as a major advantage of the general- 
ized harmonic formalism. Pretorius [ll. \3l\. building on 
Gundlach et al. [32j], introduced such a modification in 
his generalized harmonic evolution code, leading to the 
first ever simulation of a full orbit and merger of binary 
black holes. In Rcf. 33], this system was converted into 
an explicitly first-order, linearly degenerate, symmetric- 
hyperbolic form. An extensive body of mathematical lit- 
erature exists on systems of this form, and they are also 
very well suited to highly accurate multidomain pseu- 
dospectral collocation methods. While this first-order 
generalized harmonic system is perfectly acceptable for 
numerical relativity, and indeed is now used for nearly 
all simulations currently being done by the Caltech and 
Cornell numerical relativity groups, there would be con- 
siderable value in implementing a KST system of compa- 
rable stability. For one thing, the KST system involves 
just over half as many fields as the first-order generalized 
harmonic system, so it could provide a considerable im- 
provement in code runtime. The KST systems are also 
closer to the evolution systems used historically in nu- 
merical relativity. Gauge is specified again in terms of 
(densitized) lapse and shift. So a large body of research 
on gauge conditions can be more easily applied. 

As the standard KST systems already have adjustable 
parameters, it is an interesting question whether any 
of them can have constraint-damping properties. Even 
without calculations, it is clear that the answer must be 
"no." With the generalized-harmonic system presented 
in Ref. [33| . two parameters are available, 70 and 72, 
which tune the stability of small constraint-violating per- 
turbations of flat spacetime. The inverses of these pa- 
rameters (up to factors of order unity) define timescales 
on which short-wavelength constraint modes will decay 
exponentially. All of the free parameters of the stan- 
dard KST systems are dimensionless, so they cannot fix 
any preferred timescale for the damping of perturbations 
of flat spacetime. Of course in the full nonlinear phase 
space, the dimensions necessary for such a timescale can 
be provided by nontrivial features of a particular solu- 
tion, for instance the mass of a black hole. Indeed, in 
Ref. 1341 ] . it was demonstrated that careful fine-tuning 



of the parameters can considerably extend the lifetime 
of a black hole simulation. This effect cannot really be 
considered constraint damping, though, as the optimum 
choice of parameters depends strongly on the details of 
the initial data, and in fact it must fail for small fea- 
tures in asymptotic regions, where the above flat space 
arguments begin to apply. 

In the present paper, the KST systems will be gener- 



alized even further, introducing one new free parameter 
with dimensions of inverse-time, that can be considered a 
mechanism for constraint damping. This generalization 
is directly analogous to methods used in the past [33| 
for controlling the stability of the constraints that ap- 
pear when an evolution system that is second-order in 
spatial derivatives is reduced to first-order form by the 
introduction of auxiliary fields. Such a constraint exists 
in the KST systems, and while our methods are only de- 
signed to control this particular constraint, the intricate 
coupling of the constraints, in their evolution, extends the 
constraint damping effect to all (finite wavelength) con- 
straint modes, including the Hamiltonian and momentum 
constraints. 

The format of this paper is as follows. In Sec. [TT] the 
intuitive mechanism behind the new constraint damping 
terms is sketched out in the context of a simple toy model. 
In Sec. IIII1 analogous modifications are applied to the 
five-parameter KST evolution systems, and conditions on 
these modifications are noted to keep the resulting system 
hyperbolic. In Sec. HVl the hyperbolicity of the constraint 
evolution system is investigated. In Sec. [V] the effective- 
ness of these modifications on constraint-violating per- 
turbations of flat spacetime is seen in detail. In Sec. IVI1 
the available parameter freedom is summarized, and in 
Sec. IVII1 results are presented of a few simple numerical 
simulations using this new evolution system. 



II. ILLUSTRATION OF A SIMPLE MODEL 
SYSTEM 

Before jumping into the full equations of general rel- 
ativity, it would be instructive to outline the constraint 
damping idea in the context of the simplest hyperbolic 
system, the scalar wave equation: 



rf u dM = 0, 



(2.1) 



for a real scalar field tp, where rj^ is the Minkowski met- 
ric in Cartesian coordinates. This equation involves only 
one field, and no constraints, but it is second-order in 
both space and time. The second-order derivatives are 
removed by promoting the first derivatives of ip to inde- 
pendent fields. The first time derivative (up to a con- 
ventional minus sign) will be denoted with the symbol ir. 
The wave equation now becomes the system: 



d t ifj = -7T, 

<9 t 7r = -S lj didjip, 



(2.2) 
(2.3) 



where latin indices refer to the spatial coordinates of the 
chosen inertial frame. The one- form <pi, defined by a 
new constraint can be used to remove second spatial 
derivatives from the system: 



d t ip 
d t ir 
C, 



di4> -<f>i = 0. 



(2.4) 
(2.5) 
(2.6) 
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Before this can be useful for free evolution, an equation 
is needed for evolving </>j. This equation is commonly 
derived by equating the spatial coordinate derivatives of 
both sides of Eq. (|2.4[) . commuting partial derivatives on 
the left, and substituting the constraint to find: 

dtfa = -di-K. (2.7) 

This equation closes the evolution system, and leaves us 
with a nice first-order symmetric-hyperbolic representa- 
tion of the scalar wave equation. However the newly 
defined constraint field is only marginally stable. As is 
easily verified by direct substitution of the constraint def- 
inition and the evolution equations, the constraint is now 
conserved: 



system at principal order, and is clearly symmetric- 
hyperbolic. The existence of a positive-definite sym- 
metrizing inner product is independent of the basis of 
dynamical fields. Indeed, the obvious symmetrizer for 
the transformed system, 

dS 2 = Adijj 2 + dn 2 + S lj d<p l d(/) 3 , (2.15) 
= AdV 2 + (dn - 7dV) 2 + S i: > dfadfa, (2.16) 

when expressed in terms of 7r, is positive-definite (for 
positive A) and symmetrizes the untransformed system. 
Therefore, this constraint-damped form of the scalar 
wave system is symmetric-hyperbolic for any fixed choice 
of the damping timescale. 



d t Ci = 0. 



(2.8) 



This shows that exact constraint satisfaction should be 
preserved within the domain of dependence of the initial 
data, but also that any violations that may arise will be 
preserved as well. 

Because the constraint is linear in undifferentiated <f>i 1 
anything added to the right side of Eq. (|2.7|) will transfer 
directly to the evolution equation implied for the con- 
straint. For example, if the equation is changed to: 



di c 



-di-K 



■j(diip - (pi), 



(2.9) 
(2.10) 



for some constant 7, then the constraint-satisfying solu- 
tion space is unchanged, but the evolution of the con- 
straint becomes 



d t Ci = -'fCi 



(2.11) 



Thus, with this method, the constraint can be damped 
(assuming hyperbolicity is preserved) exponentially on 
an arbitrary fixed timescale 7 . As we will see when we 
discuss the KST system, the constraint damping effect 
can extend even farther than the constraint that appears 
in the reduction. Even those constraints that exist before 
the reduction to first-order form can be damped. 

The question of whether this modification preserves 
the clear symmetric hyperbolicity of the standard reduc- 
tion is important, and a very simple argument shows that 
hyperbolicity is not affected. If a linear change of vari- 
ables is made (in other words, a change of basis on the 
vector bundle of dynamical fields), defining tt := ir — jtp, 
then all modifications of the principal part of the funda- 
mental evolution system disappear: 



<9*7T 

d t (/>i 



0. 



-&7T 



(2.12) 
(2.13) 
(2.14) 



(the symbol "~" means that all nonprincipal - in this 
case, algebraic - terms have been omitted). This trans- 
formed system is exactly the same as the unmodified 



III. THE MODIFIED KST EVOLUTION 
SYSTEM 

The Kidder-Scheel-Teukolsky [25| evolution equations 
are a five-parameter 1 generalization of the standard first- 
order representation of the classic ADM [26j equations, 
in the form advocated by York [27| : 



dtgij-Crfgij = -2NK i: 



(3.1) 



dtK, 



NR tj + N{KK l3 - 2K l k K, 
-ViV ? -iV. 



kj , 



(3.2) 



The dynamical fields are {gij, Kij}, the metric intrinsic 
to the slice of constant t, and the extrinsic curvature of its 
embedding in spacetime. The gauge fields {N, N 1 }, lapse 
and shift, determine the evolution of the coordinates. 

The Ricci tensor Rij written above is that of the spa- 
tial metric g^j, so it implicitly involves second spatial 
derivatives of g+j. The evolution system can be reduced 
to first-order form by promoting the partial derivatives 
of the spatial metric functions to an independent (non- 
tensorial) three-index field: 



1 



-A 



Wij- 



(3.3) 



As long as D^ij, under its own evolution, properly rep- 
resents dkgij /2, it can be substituted for any derivatives 
of gij. This renders the ADM system first-order. 

This evolution system, like that for the scalar field, de- 
scribes physics only when certain constraint fields vanish: 



C 

Ci 



(R-K^K 11 +K 2 ) 



1 

2 

^(Kij-Kgij), 
dkgtj — 2Dkij- 



(3.4) 

(3.5) 
(3.6) 



A twelve-parameter system also exists, employing redefinitions 
of the fundamental dynamical fields and associated constraint 
substitutions. Here we will ignore this extra freedom. 
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The Hamiltonian and momentum constraints, C and Ci, 
must vanish (in vacuum) throughout each spatial slice, 
according to the four Einstein equations not represented 
in Eq. (|3.2[) . The three-index constraint Ckij vanishes 
when Dkij properly represents d k gij/2, in analogy with 
the constraint of the scalar field system. 

The new field, D k ij, is to be considered independent 
in free evolution. An evolution equation must be defined 
for this field, one that is consistent with the satisfaction 
of the three-index constraint above. The equation anal- 
ogous to Eq. ([27T]) is 

d D kij = -d k (d Q gij), (3.7) 
= -dkiNKij), (3.8) 

where the shorthand do refers to the derivative, dt — Cfi, 
along the normal to the spatial slice 2 . 

The evolution system has now been written in first- 
order form, and we can begin to ask about its hyperbol- 
icity. Kidder, Scheel and Teukolsky [25| have shown that 
the above system can be rendered strongly hyperbolic 
with a few simple modifications. The first of these is com- 
monly referred to as densitization of the lapse. Rather 
than fixing N directly, we fix a related field Q defined by 

jV = 5 T°exp(Q), (3.9) 

where g is the determinant of <?jj and 70 is a constant 
nonzero parameter. The occurences of dkgij that then 
arise arise from the V.;VjiV term in Eq. (|3.2p are then 
replaced by 2D k ij . The second modification required for 
hyperbolicity is the addition of terms to the evolution 
equations for K+j and Dhij , that are proportional to the 
constraints: 

doKij = ... + liNg l] C + l2 Ng ab C a{l])b , (3.10) 
d D kij = ... + -lzNg k (fi ]) + -j 4 Ng lj C k , (3.11) 

where the ellipses refer to the right sides of Eqs. (|3.2[) 
and (|3.8p . The four-index object C k uj ■— 2di/ l .Dn i j used 
here can be thought of as another constraint, but it van- 
ishes automatically whenever the three-index constraint 
vanishes. The parameters {70,..., 74} do not affect the 
physical solution space of the equations in any way, but 
they directly affect the principal part of the evolution 
system. In Ref. (25[, Kidder, Scheel and Teukolsky deter- 
mined sufficient conditions for these parameters that ren- 
der the evolution system strongly hyperbolic. In Ref. 
and an appendix of Ref. [3^, these arguments were ex- 
tended, and it was also made clear on what subset of 
the parameter space the equations satisfied the stronger 
condition of symmetric hyperbolicity. 



2 Note that do, involving a Lie derivative, commutes with the par- 
tial derivative d k - It is this commutation that defines the action 
of the Lie derivative on the nontensorial field D ki j . 



The focus of this paper is a further modification of 
Eq. p. lip along the same lines as that described in 
Sec. HU Here the goal is to modify the evolution of the 
three-index constraint, which in the ordinary KST sys- 
tem is implied to evolve as 

d Ckij = -73-^Sffc(iCj) - JiNgijCh- (3.12) 

Note that the need for damping here is more dire than 
in the scalar field case. Hyperbolicity requires 73 and 74 
to be nonzero, so any violation of the momentum con- 
straint will feed directly into the three-index constraint. 
This is countered as in the previous section by includ- 
ing terms proportional to the three-index constraint in 
the evolution equation for D k ij ■ As we will see in a mo- 
ment, multiples of the traces, denoted C\ := g v C k ij and 
C'j := g ht Ckij, must be added in separately, so the result- 
ing evolution equation is: 

d D klj = ... + ~K3Ngk(iCj) + ^lANgijCk 

1 2 1 , 

+ -N l7 C (l g ])k + -N-y 8 C {i gj) k 

+lN l9 C 2 k9ij . (3.13) 

The term proportional to 75 is analogous to the term 
proportional to 7 in Eq. (|2.9p . The terms proportional 
to 76, 77, 78 and 79 are necessitated by the hyperbolicity 
conditions, which we now consider. 
The principal part of this system is: 

d o9ij ~ 0, (3.14) 
BoKii c N[g ab S^ - (1 + 72)5^^ 

-(l- l2 )g bc S^8f ) + (l + 2 l0 )g c %8 b j) 
~lig ad g bc g* 3 +lig ah g cd ^\daD bcd , (3.15) 

-^74.9 ca rf + ^733 ab 5 fe (^ c ) 
+ \l49 ab 9^ c k }d c K ab 
+N[^ 5 S c k 6tS b + \^5tg ah 9ij 

1 ca cb 1 ab re 

+ 2775 0{i9j)k + 7^785 \i9j)k 

+ \l99 ca S b k9l3 }d c g ab . (3.16) 

Hyperbolicity is well established for the standard KST 
system, 75 =76 = •■■ =79 = 0, so as in Sec. UH it is 
best to seek a linear change of variables to reduce the 
constraint-damped system to the standard system. Con- 
tinuing the analogy with the scalar field, define 

Kij := K i:j - -mgij- (3.17) 
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Then the equation for doKij has the same principal part 
as that for d^Kij. The equation for doD k ij becomes 



IV. HYPERBOLICITY OF THE CONSTRAINT 
EVOLUTION 



1 



d D kij ~ -N[6° k S?6S - -j 3 g ca g k{i 6 b } 



+2-745° b g i j6l\d c K ah 
+N[\{l & -\l i l,)5lg ab g l3 



1 



1 



+ 2(77 + ^l3l 5 )g ca S b {i g j)k 
+\(ls - \izi5)g ab 5 c (% g 3)k 
+ ^(79 + \l^)g ca S b k9l] ]d c g ab . (3.18) 
If, for arbitrary 75 , the further parameters are fixed as 



76 = 


1 

7J7475 


(3.19) 


77 = 


1 

-^7375 


(3.20) 


78 = 


1 

7J7375 


(3.21) 


79 


1 

-"2-7475, 


(3.22) 



then the principal system, in the transformed variables, 
loses any reference to 75 , . . . , 79 . 

dogtj a 0, (3.23) 
doKij ~ N[g ab 5°5* - (1 + 72)5°"^^) 

-(1 - 12)9^8% + (1 + 2^ )g c %5 b j) 

-ng ad g bc gv+iig ab g cd g l j]daD bcd , (3.24) 

d D kij ~ -JV[^a§ - ^73ff co 3fe(^)) 

-7^745^4 + J 7 3<? ab <Mi) 



+ 2745° b g ij S k : }d c K ab . 



(3.25) 



This is the principal part of the standard KST system. 
So for any value of the parameter 75, with 76, ...,79 fixed 
by Eqs. (|3~T9]) - (|3~22]) 3 . the hyperbolicity of our modified 
system is the same as that of the corresponding standard 
KST system. 



3 This requirement can be weakened somewhat. There is one fur- 
ther degree of freedom, shared between 76 and 78, which will still 
preserve hyperbolicity. However when this degree of freedom is 
utilized, the simple argument used here must be replaced either 
by a somewhat more subtle argument, or a significantly more la- 
borious one. We have not yet found a use for this further degree 
of freedom, so here we restrict attention to the simpler case. 



Let us now turn our attention to the evolution of 
the constraint fields in our modified KST evolution sys- 
tem. Given the definitions of the constraints in terms of 
the fundamental dynamical fields, an evolution system 
for the constraints follows from our fundamental evolu- 
tion equations. It is important, for the construction of 
constraint-preserving boundary conditions, that this sys- 
tem also be symmetric- hyperbolic |30L l35l . l36| . The prin- 
cipal part of this evolution system can be expressed as, 

d a C ~ -i(2 -73 + 274)Ay j d J C J 
+ri(78 -276-75W^d 



+-(77-279+75)AWq 
-(1 + 2 7l )iV0 i C 



(4.1) 



+ -Ng kl g ab [(l - l2 )d k C laU + (1 + l2 )d k C, 



-(l + 2f )d k C Hab \, 



doC k ij ~ 0, 

doCabij — — 



ailb 

(4.2) 
(4.3) 



-^Nj 3 g l[a d b] Cj - -Nj 3 g j[a d b] Ci 
-Nj4gijd[ b C a ] - N^gijdihC^ 
^N-y 7 g l[a d b] C^ - ^N-f 7 g j[a d b] Cf 

-^N 1% g l[a d b] C] - -Nj 8 g j[a d b] el 
-Nl9gijd [b C 2 a] , 



(4.4) 



where the four-index object C k uj :— 2dy k Di^ is consid- 
ered an independent constraint, so that the constraint 
evolution system is first-order. 

The inclusion of the terms proportional to 75, ...,79 
in Eqs. I|4.ip - (|4.4[) has seriously complicated this sys- 
tem. Let us consider, however, the case considered above, 
where 75 is arbitrary, and the other parameters are fixed 
by Eqs. (j3~T5)) - (j3~2"2"j) . In this number of remark- 

able simplifications occur and the above system can be 
written as 



d C 
d Ci 



doC k ij 
doC ab ij 



--(2- 73 + 2 7i )Ng»d i C j , (4.5) 
-(1 + 2j 1 )Nd i C 

+\Ng kl g ab [{l - li)d k C laU + (1 + j 2 )d k C ailb 
-(l + 2 l0 )d k C ltab ], 

0, 

-2 N 73gi[ a db]Cj ~ 2 Nl39 rt adb & 



-N^g^dihCa] 



(4.6) 
(4.7) 

(4.8) 



6 



defining the new combination Ck ■= C k + 575^ — 575^. 

This constraint evolution system has the same princi- 
pal part as the standard KST constraint system. Thus, 
when the parameters are chosen by Eqs. (|3.19[) - (|3.22[) . 
the hyperbolicity of the fundamental and constraint evo- 
lution systems are independent of the parameter 75, so 
our modifications do not alter the hyperbolicity of these 
systems. 

V. STABILITY OF CONSTRAINT FIELDS 
UNDER FREE EVOLUTION 

Analyzing the stability of the constraint evolution sys- 
tem in generic simulations is essentially no different than 
the full numerical relativity problem itself. In order 
to get some handle, at the analytical level, on the ef- 
fect of our modifications, we consider constraint violat- 
ing perturbations of Minkowski spacetime. Obviously 
these estimates will not be completely relevant in sim- 
ulations of interest, but at least in the limit of short- 
wavelength perturbations, the dependence on the space- 
time background should be minimal. In this sense, sta- 
bility of short-wavelength constraint-violating perturba- 
tions of Minkowski spacetime is a necessary condition for 
constraint damping in general. And while our analysis of 
long-wavelength modes may not be directly relevant for 
evolutions of curved spacetime, unstable long-wavelength 
modes should at least be disconcerting, as a signal that 
instabilities are likely in general simulations. 

This analysis involves the full (not just principal) con- 
straint evolution system, linearized about the limit that 
gij = %, Kij = D klj = 0, N = 1, N* = 0. In this 
context, the full constraint evolution system becomes: 

d t C = ~(2- 73 +2 74 )5« [diCj + ^diCj 

(5.1) 

dtC, = -(l + 2 7 i)^C 

+l S kl 6 ab [(l - l2 )d k d [a C l]bl + (1 + l2 )d k d Vl C a]lb 
-(l + 2 7 o)d fe %e /]afc ], (5.2) 

dtCkij = —j5Ckij 

— 73<5fc(iCj) - 74%Cfc 

-^73 7 5<**(<C]) - ipAlhhfil 

+ 2'Y3l55 k (iCj) + 7;7475%Cfc- (5.3) 

Notice that we are no longer considering Ckiij an inde- 
pendent constraint field. In actual evolutions, where the 
fundamental fields are evolved, not the constraints, the 
three- and four-index constraints satisfy the identity 

Ckiij = d[iCk]i r (5.4) 
Violations of this identity will not appear in evolutions. 



Now the above system is simplified by resolving all 
constraint fields into Fourier modes. This has the formal 
effect of replacing all spatial derivatives dj with — ikj, an 
imaginary unit times a propagation vector kj . The result 
is a system of coupled ODEs for the various constraint 
modes c A (ki, t): 

8 t c A = M^c B . (5.5) 

Each eigenvector of M^{ki) evolves as exp(si) for some 
s(fcj). The real part of s is the rate of exponential growth 
(or damping, if negative) for the corresponding mode. 
Due to the rotational invariance of the problem, these 
eigenvalues should depend only on the magnitude of fej, 
so the propagation vector is decomposed as ki — kvii 
where m is a unit vector. 

This eigenvalue problem naturally reduces into sub- 
spaces according to various possible spin weights about 
the propagation direction rij. There is a five-dimensional 
space of longitudinal modes: {C, C n , C„, C„, C nnn }, where 
C n := n l Ci, etc. There is also a five-dimensional space of 
transverse vector modes: {Ci,C},C],Ci nn ,C nn i}, where 
capital latin indices now refer to a two-dimensional vector 
basis orthogonal to n l . The remaining constraint fields, 
with higher spin weight, are represented among the var- 
ious projections of the totally tracefree part of Ckij ■ A 
glance at Eq. (|5.3p shows that all of these high spin- 
weight fields propagate trivially with s = —75 indepen- 
dent of wavelength. They are therefore damped expo- 
nentially on the timescale 7^ for positive 75. The lon- 
gitudinal and transverse constraint modes require more 
careful consideration. 



A. Transverse vector constraint modes 

The growth rates of the transverse vector modes are 
related to the eigenvalues of a five-by-five matrix. Three 
of these eigenvalues simply equal —75. The remaining 
two are solutions of a quadratic equation, and depend on 
wavelength as: 

s{k) = ~ l5 T±^j\^-vlk\ (5.6) 
where we define the convenient shorthand 

r:= i(2 - 73 + 2 74 ), (5.7) 

and V2 is one of the characteristic speeds of the KST 
system (relative to hypersurface-normal observers), 

v i : = ^73(1 - 372 - 470) - ^74(1 + 670). (5.8) 

Notice that one mode is undamped in the long- 
wavelength (k — > 0) limit, where one root in Eq. (|5.6|) 
becomes zero. This is not surprising: other constraint- 
damped representations of the Einstein system have the 
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same property [23j, [32], l33[. In practice, long wave- 
length constraint modes should be killed off by proper 
constraint-preserving boundary conditions. 

In the short wavelength (k — > oo) limit, the dispersion 
relation becomes 



s(k) — > -^75r ± iv 2 k. 



(5.9) 



These represent propagating modes of the constraint 
system, damped at short wavelength on the timescale 
(575]?) -1 . Notice the significance of the constant V. 
Most of the modes require 75 > for damping, so the 
damping of the modes referred to in Eq. (|5.9j) requires 
r > as well. Thus, the damping condition places a new 
restriction on the standard KST parameters {70, ...,74}, 
beyond the conditions they must satisfy for the system 
to be hyperbolic. 



B. Longitudinal constraint modes 

The longitudinal modes again involve the eigenvalues 
of a five-by-five matrix. In this case two of the eigenval- 
ues are simply —75. The rest are the roots of the cubic 
polynomial 

s 3 + 75IV 2 + k 2 vls + /c 2 75 r(l + 2 7 i) = 0, (5.10) 
where 1*3 is another characteristic speed, given by 



«§ ■=l( 1 + 2 7i)(2 - 73 + 2 74 ) - ^7273- (5.11) 



Rather than giving complicated analytic expressions for 
the roots of this polynomial, we simply consider asymp- 
totic limits in k. First, in the long-wavelength (k = 0) 
limit, two roots vanish and the third is — 7 5l\ This is 
very similar to the long-wavelength behavior of the vec- 
tor modes. 

In the short-wavelength limit, the polynomial becomes 
singular. The terms proportional to k 2 dominate the 
polynomial, leaving a linear equation. The root of this 
linear equation, s = — i>3~ 2 75r(l + 2 7 i), is the regular root 
of the polynomial in this limit. The two remaining roots 
disappear in the limit k — > 00. These singular roots cor- 
respond to traveling modes, with imaginary part linear 
in k in this limit. They can be found by substituting for 
s a power series in fc, s — s±k + sq + + ... in the 

above polynomial and solving the resulting polynomial 
order-by-order in k for the coefficients Si. The result is 

s{k) = - ^T 21 ) ± "0* + Oik- 1 ). (5.12) 

So the damping of the traveling longitudinal modes re- 
quires that v 2 > (1 + 271) when the transverse modes are 
damped as well. 

In summary, the damping of short-wavelength 
constraint- violating modes requires that the rates 



T-2 



1 

2 75F 

U3 - 2 7 5 r(i + 271) 
1 

2 



^3~ 2 75l> 3 2 " 1 " 2 7 l) 



(5.14) 
(5.15) 
(5.16) 



be positive, where T is defined by Eq. (|5.7[) and U3 by 
Eq. (IBTTTIl . 



VI. CHOOSING PARAMETERS 

Before proceeding with numerical tests, values must be 
fixed for the free parameters. The parameters associated 
with the constraint damping terms are reasonably well 
set. The overall damping timescale is set by I/75, and 
this can be chosen to be any positive number. The other 
new parameters are determined by Eqs. (|3.19[) - (|3.22[) . 
The original KST parameters should be chosen in accord 
with hyperbolicity conditions for the fundamental and 
constraint evolution systems, as well as the conditions 
that the damping rates of Eqs. (|5.13[) - (|5.16[) be positive. 

The hyperbolicity conditions are quite complicated 
when considered in full generality To make the situation 
more tractable, here we restrict attention to the subset 
of parameter space in which all characteristic speeds are 
equal to zero or unity, relative to hypersurface-normal 
observers. The hyperbolicity conditions in this subset of 
the parameter space are spelled out in Appendix B of 
Ref. [3(|, following work in Ref. [29| . The parameters 
70, 73 and 74 are fixed in terms of 71 and 72 by the the 
conditions on the characteristic speeds: 



7o = 



73 



74 



4 72 + (5 + 3 72 )(l + 27i) 

I-72- (l + 27l)(5 + 372) 

472 + (5 + 37 2 )(l + 2 7l ) 



(6.1) 
(6.2) 

(6.3) 



The fundamental evolution system is then symmetric- 
hyperbolic so long as the following inequalities are satis- 
fied: 



f < 72 < 



4 72 + (1 + 2 7l )(5 + 3 72 ) ^ 0. 



(6.4) 
(6.5) 



Constraint damping requires that the rates rj of 
Eqs. (|5.13[) - (|5.16p be positive. This in turn requires 
that T > and that 



< 1 + 2 7 i < vi 



(6.6) 



For the numerical simulations presented in the next sec- 
tion, 71 = -1/4. The parameter T := (1/2) (2 -73 + 274) 
can be expressed in terms of 72 using the above expres- 
sions for 73 and 74: 



ro := 75 



(5.13) 



5 + 3 7 2 



10 + 672 

4 7 2 + (5 + 372)(l + 2 7l ) ~ 5 + ll7 2 ' 



(6.7) 
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where the last equality is restricted to the case 71 = 
— 1/4. In the allowable region for 72, T can be set equal 
to any value greater than 2. Here we choose T = 5/2. 

The various parameters, and the associated growth 
rates, come out to: 



7o 


1 

2 


(6.8) 


71 


1 

~4' 


(6.9) 


/ ^ 


5 

^43' 


(6.10) 


73 


43 

~ ~10' 


(6.11) 


74 


52 

~80' 


(6.12) 




= 75, 


(6.13) 


ri 


5 

= ^75, 


(6.14) 




5 

= ^75, 


(6.15) 




5 

= o75- 


(6.16) 



These parameters satisfy all the necessary conditions 
for constraint damping in perturbations of flat space- 
time, as well as those for symmetric-hyperbolic propaga- 
tion of the fundamental evolution fields. Unfortunately, 
these parameters do not satisfy all of the necessary con- 
ditions for symmetric-hyperbolic constraint propagation. 
In Ref. [iOl it was shown that when the adjustable char- 
acteristic speeds are all set to unity, the symmetric hy- 
perbolicity conditions on the fundamental and constraint 
evolution systems collude to require that 1 + 271 < 0, a 
direct conflict with our damping conditions. Unfortu- 
nately, this conflict does not appear to be an artefact 
of our condition that all adjustable characteristic speeds 
are equal to one. Monte Carlo searches over the en- 
tire available parameter space have not provided us with 
any examples of systems with constraint damping along 
with symmetric- hyperbolic propagation of the fundamen- 
tal and constraint fields. 

In principle, this conflict is very serious. At timelike 
boundaries of the simulation domain, conditions must be 
imposed on fields entering the computational grid. These 
boundary conditions should be compatible with the con- 
straint equations. In Ref. [30(, such boundary conditions 
were presented. These conditions control the growth of 
a certain norm of the constraint fields. In the case of the 
parameters used here, this norm is not positive-definite, 
so control of the norm does not necessarily imply control 
of the constraint fields themselves. 

In practice, the damage done by this conflict can 
only be assessed with numerical simulations. While 
the constraint evolution is not symmetric-hyperbolic, it 
is strongly hyperbolic, so the boundary conditions of 
Ref. [30( can still be applied, even if they may not have all 
of the desired effects. In fact, the numerical results of the 



following section demonstrate that constraint-preserving 
boundary conditions are quite effective in these simula- 
tions. Perhaps this can be explained heuristically by the 
fact that the "timelike" degree of freedom in the con- 
straint evolution (the one whose violations could com- 
pensate, in the indefinite norm, for violations of the 
other constraints) is very well controlled by the constraint 
damping. 

It should also be noted that without the constraint 
damping terms, the particular parameter set used here 
leads to very unstable evolutions. In the following sec- 
tion, we will not make comparisons with the undamped 
case, 75 = 0, as those cases immediately become unsta- 
ble. This could be due, in part, to the lack of symmetric- 
hyperbolic constraint evolution. At any rate, when the 
constraint damping terms are included, the evolutions 
become remarkably stable. 



VII. NUMERICAL TESTS 

The following numerical tests were carried out using 
the Spectral Einstein Code developed over the last few 
years by the numerical relativity groups at Cornell and 
Caltech. The code uses multidomain pseudospectral col- 
location methods to resolve the fields in space with ex- 
ponential accuracy. Integration in time is implemented 
by the method of lines, using in this case a fourth-order 
Runge-Kutta scheme. More details on this code and its 
remarkable accuracy can be found in Ref. [37j and refer- 
ences therein. 

The spectral representation of the computed fields 
is done in accordance with the topology of the spa- 
tial domain. The present simulations are of a single 
Schwarzschild black hole, in Kerr-Schild [38| coordinates. 
The spatial domain is made up of a family of concentric 
thick spherical shells. The fields are therefore resolved 
into spherical harmonics in the angular directions, mul- 
tiplied by Chebyshev polynomials in the radial direction. 
The innermost boundary is inside the black hole hori- 
zon, so no boundary condition is needed there. At the 
outermost boundary, the constraint-preserving boun dary 
conditions presented in Ref. [30] are used. As in Ref. [30] , 
tensor spherical harmonic components of the four highest 
I values are discarded after each time step. No filtering 
appears to be necessary in the radial direction. 

Figures [1] and [5] demonstrate the stability and expo- 
nential convergence of these simulations. Figure [I] is a 
plot of the error norm: 

I H| 2 := J (Sg^Sgij+SK^SKij+dD^SDk^dV, 

measuring the difference between the computed solution 
and the reference Kerr-Schild geometry. Figure [5] shows 



9 



a positive-definite norm of the constraint fields: 

||C|| 2 := J [c 2 + ^C% + ^C kij C kij + ^C mj C m ^jdV. 

(7.2) 

We normalize these quantities, dividing by norms that 
involve similar terms, but that should not be expected to 
vanish. The error norm is divided by the overall solution 
norm: 

I Ml 2 : J (/y''//.., + A "A.., + D ki W kij ) dV, (7.3) 

and the constraint norm is divided by a similar norm of 
the first derivatives of the computed fields: 

\\du\\ 2 := J{g kc 9 ia 9 3h d k9l3 d c g ab 
+g kc g ia gi h d k K lJ d c K ab 

+g ld g kc g ta g 3b diD ktJ d d D cab )dv. (7.4) 

All indices are raised and lowered with the computed 
metric gij. 

Here, the inner (excision) boundary is at 1.9M, and the 
outer boundary is at 41. 9M. This domain is divided into 
eight subdomains, each of coordinate thickness 5M. This 
is the same domain used in Ref. (jjjj. Note that the con- 
vergence stops at the highest resolution presented here, 
overtaken by exponential growth that is not yet ap par ent 
in the constraint fields shown in Fig. [2] In Ref. [30| , a 
"gauge instability" was mentioned, associated with one 
particular boundary condition. Presumably, this is the 
same instability apparent in Fig. [I] in which case it could 
be expected that convergence would improve as the lo- 
cation of the outer boundary is moved farther into the 
asymptotic regime. 

As a test of this hypothesis, the highest-resolution run 
in Fig. Q] was repeated on larger domains, keeping resolu- 
tion fixed but adding extra subdomains to place the outer 
boundary at coordinate radii 61.9 M, 81.9 M, 101.9 M. 
Figure [3] demonstrates the improvement in the overall er- 
ror norm. Least-squares fitting of the data in that plot 
show that the late-term growth in this error occurs ex- 
ponentially on a timescale proportional to the square of 
the coordinate position of the outer boundary. Figure 0] 
shows the growth of constraint energy in these simula- 
tions. Until an exponential instability sets in, apparently 
triggered by the overall loss of accuracy of the simula- 
tion, the constraint fields grow roughly as the square 
root of coordinate time. On the largest domain, this 
slow growth persists beyond 15,000 M, when exponen- 
tial growth takes over at a rate that would allow the 
simulation to survive until nearly 50,000 M. 

It is of some interest to verify the effectiveness of 
the constraint-preserving boundary conditions used in 
these simulations. As noted in the previous section, 
since the characteristic matrices of the constraint 
evolution system are symmetric only with respect to 
a Lorentzian norm, there is no reason to expect these 




l/M 



FIG. 1: Norm of the error ||(5u||/||it||, relative to the reference 
solution, on a fixed domain extending from minimum coordi- 
nate radius 1.9M to maximum 41. 9M. The domain is broken 
into eight shells each of thickness 5M and radial resolution 
N r , chosen on four different runs as N r = 8, 11, 14, 17. The 
constraint damping terms presented in all of these simulations 
have 7s = 0.6A/ -1 . 
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FIG. 2: Constraint norm ||C||/||9u||, for the same runs plotted 
in Fig. [TJ 



conditions to control the influx of constraint violations. 
In Fig. El the simulation with R ma x = 41.9 and N r = 17 
(in each subdomain) is repeated using conventional 
boundary conditions. These boundary conditions freeze 
the incoming characteristic fields of the fundamental 
evolution system to their initial values. These "freezing" 
boundary conditions control a positive-definite norm of 
the fundamental evolution fields, so the initial boundary 
value problem is known to be well-posed by standard 
theorems. However, the figure clearly demonstrates the 
superiority of the constraint-preserving boundary condi- 
tions in this context, not only for constraint satisfaction, 
but for overall stability. Perhaps the effectiveness of the 
constraint preserving conditions is not robust, perhaps 
it will fail when the conditions are applied in more 
dynamical spacetimes. This possibility is an important 
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1/M 



FIG. 3: Error norm ||<5u||/||w|| for runs in the constraint- 
damped system with outer boundary at r = 41.9 M, 61.9 M, 
81.9 M, 101.9 M. The long-term growth of the error norm 
occurs exponentially on a timescale proportional to the square 
of the coordinate position of the outer boundary. 
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FIG. 4: Constraint norm |C||/||<9u|| of the same runs as those 
in Fig. [3] The constraints grow as t 1 ^ 2 until eventually driven 
exponentially by the overall loss of accuracy demonstrated in 
Fig. E 

avenue for further investigation - if this stability is found 
not to be robust, then either the spatial domain will 
need to be compactified to remove timelike boundaries, 
or further modification of the KST system will be needed 
to combine the constraint damping effects outlined here 
with truly symmetric hyperbolic constraint propagation. 



VIII. DISCUSSION 

A generalization of the five-parameter KST systems 
was introduced, for use in numerical relativity. The 
added parameter 75 supplies a timescale on which ex- 



ponential damping can occur (or growth, if parameters 
are not chosen carefully). The hyperbolicity of the funda- 
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FIG. 5: Constraint norm ||C||/||9u|| of a black hole simulation 
with Rmax = 41.9, N r = 17 (in each of eight subdomains), 
and two different boundary conditions. 



mental and constraint evolution systems is not changed 
by this modification, but the effect that the constraint 
damping has on perturbations of flat spacetime is partly 
dependent on the same parameters that determine hy- 
perbolicity. Parameters can be chosen such that all con- 
straint modes are stable in perturbations of flat space- 
time, but not when the constraint fields are required to 
evolve in a symmetric-hyperbolic manner. Nevertheless, 
single black hole simulations using constraint-preserving 
boundary conditions are convergent, and what instabili- 
ties exist appear to be dominated by constraint-satisfying 
modes, associated with a gauge instability in the outer 
boundary condition. 
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